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^ ! Abstract 

. Correlations between the arrival time and the energy of photons emitted in outbursts 

^ J> \ of astrophysical objects are predicted in quantum and classical gravity scenarios and 

can appear as well as a result of complex acceleration mechanisms responsible for 
the photon emission at the source. This paper presents a robust method to study 
such correlations that overcomes some limitations encountered in previous analysis, 
and is based on a Likelihood function built from the physical picture assumed for the 
emission, propagation and detection of the photons. The results of the application 
of this method to a flare of Markarian 501 observed by the MAGIC telescope are 

\ presented. The method is also applied to a simulated dataset based on the flare 

of PKS 2155-304 recorded by the H.E.S.S. observatory to proof its applicability to 
complex photon arrival time distributions. 
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1 Introduction 



Observations of the electromagnetic radiation from astrophysical sources ex- 
hibiting fast flux variations -mainly active galactic nuclei (AGNs) and gamma- 
ray bursts (GRBs) and pulsars- make possible the study of correlations be- 
tween the energy and the arrival time of individual photons. Two different 
effects can explain the appearance of such correlations: either they are source- 
intrinsic or have originated during the propagation of the radiation from the 
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source to the Earth. On the one hand, the emission mechanisms taking place 
at the source can cause delays between photons of different energies to ap- 
pear. In the case of pulsars, a difference in the emission location for different 
energies can cause the shape and peak time of the pulses to change with en- 
ergy [rf2] . In AGNs the delays can be caused, for instance, by the acceleration 
of the electrons responsible for the emission [3f4"] . On the other hand, if the 
emission is assumed to be simultaneous at the source for all the photons in a 
certain energy range, an energy-dependent propagation effect could also ex- 
plain the appearance of these correlations. Quantum Gravity (QG) [5foir?f8f9"] 
and classical gravitation -through the existence of wormholes PU]~ are two 
frameworks where this kind of propagation effects may show up. 

It is widely speculated that space-time is a dynamical medium, subject to 
quantum-gravitational effects that cause space-time to fluctuate on the Planck 
distance scales. Model realizations of QG predict that quantum fluctuations in 
the space-time metric make it appear 'foamy' on very short distance scales (see 
|11|I12|I13|[T4"] for reviews). The propagation of photons through this fluctuating 
space-time might induce an energy dependence of the speed of light, resulting 
in an observable difference on the propagation time for photons of different 
energies. The anomaly induced by QG will always be a small perturbation 
of the assumed light speed c, and can therefore be treated as a correction 
of the form Ac/c = —E/Mqgi or —E 2 /Mq G2 where E is the photon energy, 
Mqgi and Mqg2 are respectively the effective QG scales in the linear and the 
quadratic term, and the minus sign indicates that most of the models expect 
this effect to be subluminal. 

The possibility to observe QG efects in the propagation of photons from astro- 
physical sources was proposed by Amelino-Camelia et al. [15J. If two photons 
of different energy are emitted simultaneously from an astrophysical object, 
the expected delay between their arrival times when they are detected in- 
creases with the distance to the source and with the energy of the photons. 
Therefore, the maximum sensitivity to energy dependent propagation effects 
is expected from observations of very fast flux variations coming from sources 
at large distances that emit photons up to very high energies. 

The first candidates considered for such observations were gamma-ray bursts 
(GRB) [15] . but the first experimental result on possible energy-dependent 
speed of light came from the measurement of a flare of the active galactic 
nucleus (AGN) Mkn 421 at TeV energies by the Whipple gamma-ray tele- 
scope [in], claiming a lower limit to the energy scale of quantum gravity of 
Mqgi > 6 x 10 16 GeV . Other bounds have been obtained by studying the 
emission of pulsars [IT], and a combined analysis of many GRBs yielded to a 
robust lower limit of M QG i > 0.9 x 10 16 GeV [H]. 

The highest sensitivity of the new generation of ground-based gamma-ray 



2 



telescopes was expected to have a clear impact in these studies [19J and, in- 
deed, has provided new observations of fast flares of AGNs at TeV energies 
with richer photon statistics than the previous generation of instruments could 
[2"0~f2"T] , enhancing the sensitivity of these measurements to any energy depen- 
dent propagation effect. 

This paper presents a new method of analysis that can be applied to any 
observed set of photons to search for any kind of correlations between their 
arrival time and energy. This analysis technique is described in section [2] and 
applied to the Mkn 501 flare observed by the MAGIC telescope in July 2005 
[2D] in section [3J An application to a simulated set of photons based on the 
more complex flare structure of the PKS 2155-304 outburst observed by the 
H.E.S.S. collaboration in July 2006 [21] is discussed in section [U Finally, the 
conclusions of this work are summarized in section [5j 



2 Description of the method 

The pioneering search for QG efects using TeV gamma-ray data by the Whip- 
ple Collaboration [TB] was based on few tens of gamma-ray events that where 
both binned in time and energy. 

Given the fact that the statistics for sources at the largest distances from us 
and photons of the highest energies (conditions that maximize the expected 
photon delay) is usually scarce, in order to make an optimal use of the infor- 
mation of the arrival time and estimated energy of each recorded photon the 
analysis method used should be unbinned. 

Some analysis methods, such as the Modified Cross Correlation Function 
(MCCF) [22123] or the Continuous Wavelet Transform (CWT) [23] have been 
widely used for similar analysis [T5|25ll26f27f28j . These methods have been 
applied using the information of the arrival time for each individual photon 
without binning, but binning the energy in bands in order to determine the 
time lag as a function of the energy and therefore did not make optimal use 
of the photon information. 

In addition, as mentioned before, measurements at the highest photon energies 
maximize the expected delays between photons, but all the sources detected at 
the GeV-TeV energy range show steeply falling energy spectra, reducing the 
photon statistics as the energy increases. Some of the above analysis methods, 
as for instance the wavelet approach, [25j simply cannot deal properly with 
the low photon statistics obtained in most of the very high energy gamma-ray 
measurements. 
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Completely unbinned approaches to search for correlations in similar data 
sets are nevertheless possible. For instance the Energy Cost Function (ECF) 
method described in [2S] has been successfully used to measure a tiny time 
lag in the same data set which will be used here as an example. Nevertheless, 
that method requires the identification of well-isolated flares and cannot be 
applied directly to complex lightcurves where several overlapping flares, or 
simply other kind of features such as edges of non-contained flares, appear in 
the lightcurves. 

This work describes a completely unbinned method which overcomes the above 
limitations and can be applied in a general manner to any data set, regardless 
on the photon statistics and the lightcurve shape and structure. Moreover, it 
makes use of an optimal estimator which, if properly applied, would use all 
the information contained in the data and therefore provide with the most 
accurate estimate of the correlation between energy and arrival time. This 
estimator is the Likelihood Function built from the combined probability of 
having observed a set of photons with individual energy Ei arriving at time t%, 
and has to be constructed from a physical description of the assumed processes 
involved in the emission, propagation and detection of the photons. 

Any hypotheses about the physics mechanisms that finally produce the mea- 
sured set of photons can be accommodated in this analysis provided a reason- 
ably simple mathematical description is at hand. A possible general physics 
picture of the process is the following: 

• Photons were produced at the source following an intrinsic lightcurve 1 
F a (t s ) where t s is the time at the source and an intrinsic energy spectrum 
T(E a ) where E s is the energy of the photons at the source. It might be that 
already at the source the lighcurve depends on energy or the energy spec- 
trum depends on time. If models describing these dependencies are at hand 
for the observed astrophysical source, they can be readily incorporated in 
the approach described here. 

• These photons propagate in space and there, if there is an energy-dependent 
refraction index related to some effective energy scale parameter Mqg„, a 
delay D(E S , Mqgu, z) should affect their propagation time. 

• These photons reach the Earth and are recorded by satellite or ground- 
based instruments as photons with measured energy E and arrival time t. 
These detectors have typically an extremely good time resolution compared 
with the flare time scales but a limited energy resolution G(E — E s , cfe(E s )) 
which may be any complex function of the observed energy. 



The term "lightcurve" is widely used in astrophysics to denote the photon emis- 
sion time distribution function. 
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2.1 Probability Density Function 



A mathematical expression casting the probability density function (p.d.f.) 
describing the physics picture given above is the following: 



dP f°° 

j^j t =N J o T(E S ) C(E S , t) G(E - E s , a E (E s )) F s (t - D(E 8 , M QGn , z)) dE s 

(1) 

where 

• P is the probability density function (p.d.f.) for a photon to have observed 
energy E and arrival time t. 

• T(E S ) is the photon energy distribution at the source. For instance, for Very 
High Energy gamma rays, in most cases the energy spectrum is well de- 
scribed in first approximation by a pure power law T(E S ) = Tq ■ (E s /Eo)~ a 
where a is the differential spectral index. Nevertheless, more complex func- 
tions such as broken or curved power laws, spectral index for the photons 
coming from the flare activity different from the spectral index for the pho- 
tons coming from the steady activity, or synchrotron or inverse Compton 
spectra can be readily used. The spectrum can be obtained from a fit to 
the whole photon dataset assuming that it does not change during the flare. 
Otherwise, if the photon distribution is observed to change with time it 
can be fit to the data assuming a specific mechanism for correlation at the 
source. 

• C(E s ,t) accounts for changes in the effective area of the detector during 
the observation of the flare. In most cases this factor will be constant and 
can be neglected. It has to be included, for example, if the data comes from 
ground-based telescopes with changing observation conditions (significant 
changes in the zenith angle of the pointing, atmospheric conditions, etc.) or 
from satellites observing in survey mode. 

• G(E — E s , <7e(E s )) is the gamma energy smearing produced by the instru- 
ment. In first approximation it could be a gaussian distribution with width 
o~e{E s ) although more complex functions can also be used. This function 
has to be obtained from the study of the energy response of the detector. 

• D(E S , Mqgui z) is the eventual propagation delay as a function of the photon 
energy E s , an effective energy scale Mq Gn and the source redshift z. This 
function may accommodate linear, quadratic or any arbitrary realization 
of any possible propagation delay mechanisms. A general implementation 
of this function in terms of an effective expansion in the case of Quantum 
Gravity scenarios is discussed below. 

• t s — t — D(E S , Mqcu, z) is the actual photon production time at the source 
obtained from the measured arrival time corrected by the delay generated 
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by non-trivial propagation effects. 

• F s (t s ) is the emission time distribution of the photons at the source. For 
most astrophysical sources there are no compelling models which can be 
used to predict and/or fit the flaring lightcurves. Moreover, since the prob- 
ability P depends only on the relative delay of the photons for different 
energies, the best estimator for the time structure at the source can be 
taken to be the one of the observed gamma emission time distribution for 
the lowest energy gamma rays. In the specific case of Very High Energy 
gamma rays, this estimator is a very good choice because, since the energy 
spectra are typically steep power laws, most of the recorded gammas are 
collected at the lowest energies, and therefore the gamma statistics is then 
sufficient to infer the lightcurve as explained above. 

• N is the p.d.f. normalization. This normalization is extremely important 
for the proper use of the probability in a likelihood function. To have an 
unbiased estimation of the fitting parameters it has to be computed in such 
a way that the p.d.f. integral is equal to 1 no matter the values of the free 
parameters in the p.d.f. 

If the propagation delays are inspired in Quantum Gravity models, in a power 
n realization of Quantum gravity effects (n — 1 linear and n = 2 quadratic) 
for z << 1 we have that the propagation delay of a photon of energy E s with 
respect to the arrival time of a photon of energy E s <q is: 



E 1 } — E 7 } 



D(E s ,M QGn ,z) = ^—^-?- (2) 

M QGn H 

where z is the source redshift and H is the local Hubble constant (in inverse 
seconds). For sources at larger distances the redshift in the photon energies 
and the expansion effect in the photon path have to be taken into account, 
and then the correct expressions are: 



T1(F F y\ 1 E s-Es,Q f Z (1 + Z)dz 

D(E s ,E QG ,z) = — — / — — — (3) 



for the linear case and 



H M QGl Jo h(z) 



D(F F r\ 1 g " E -° [ Z ^ l + Z)2dz (A) 
V(E s ,Eq G ,z) = — —2 / t--t (4) 

H M^ G2 Jo h(z) 

for the quadratic case where 

h(z) = Jn A + n M {i + zy (5) 
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being £7a and VLm the standard cosmological parameters. 



2.2 Likelihood Function and fitting procedure 

With the above described p.d.f. the likelihood function for the observation of 
all the N 7 photons recorded during the flare can be build as 



In the above expression, several free parameters may appear. One of them is 
obviously the effective scale Mqgu at which propagation delays appear. Other 
parameters may try to fit the intrinsic lightcurve shape to some analytic or 
semi-analytic parametrization or the intrinsic energy spectrum. The likelihood 
function can be fitted to the data using a numerical procedure to find the set 
of free parameters that minimize the function —2log(L). 

From the conceptual point of view, when the Mq Gu parameter is changed, 
the effect is basically sliding in time the lightcurve profile as a function of 
energy to fit the possible correlation between energy and arrival time in the 
data. If there is no correlation between energy and time, or the correlation is 
not significative that would correspond to Mqcn — * ±00 and therefore Mqcn 
is not a good fitting parameter for the expected small correlations, because 
the likelihood function would be highly non-parabolic. Instead, a well-behaved 
fitting parameter is M P /MQ Gn , where M P = 2.4 x 10 18 GeV is the reduced 
Planck mass, which can have positive or negative or zero values and will likely 
show a nice parabolic behavior. Moreover, in most QG scenarios its order of 
magnitude is expected to be around 1 which makes numerical calculations 
easier. 

In the physics picture assumed, when sliding in time the lightcurve all gammas 
are expected to lead to the same "intrinsic" lightcurve, therefore all gammas, 
regardless on their energy, contribute in the fitting procedure to either deter- 
mine or test the assumed lightcurve shape even if just the lower energy ones 
were initially used to estimate the shape of the lightcurve. 

If the fit is well behaved, the fit results for all free parameters shall show a 
parabolic behavior of the log-likelihood function around the "best fit mini- 
mum". Around that minimum, an increase of A(— 2log(L)) = 1 will then give 
the single-parameter one-a uncertainties. 

Since the Likelihood p.d.f. by construction has dimensions [i?] -1 ^] -1 , its value 
at the minimum does not have a direct interpretation and does not provide a 



iV~ 

L = U 



dP 



(6) 
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simple means to estimate the goodness-of-fit for the assumed function. There- 
fore additional tests must be performed to quantify the adequacy of the physics 
picture casted in the Likelihood function to describe the data. 

One simple possibility that we have used is the following: once the Likelihood 
fit is done, the data and the Likelihood function computed at the best fit 
parameter values are binned in energy and time. Then the \ 2 value of the 
comparison of the resulting bi-dimensional histograms (in energy and time) is 
computed and compared with the number of degrees of freedom, providing a 
goodness-of-fit probability. This approach allows also an easy visualization on 
the adequacy of the fitting function to describe the data. Other possibilities 
such as integrating the likelihood function would require the use of complex 
and lengthy multi-dimensional integrations. 



3 Example of application to a Markarian 501 flare 

The blazar Markarian 501, located at z = 0.034, showed a very high activity 
in the gamma-ray band in June and July 2006. On July 9 a very strong flare in 
the TeV regime was observed by the MAGIC collaboration [20] • The analysis 
procedure explained in the previous section was applied to this set of photons 
in reference |29j together with a method based on maximizing the total energy 
in the most active part of the flare. 

In that data sample, the lightcurve was tested to be well described by a simple 
gaussian flare on top of a baseline constant in time within the observation win- 
dow (see Fig. [T]). It was tested that other mathematical functions fitting also 
well the data flare shape, such as the halving-doubling time formula used in 
the original publication [20J or even a simple triangle function, were producing 
similar final results. Therefore, the time distribution was parameterized as a 
gaussian flare of width tw and position to relative to the first gamma arrival 
time, on top of a flat background time distribution (see figured]). 

The intrinsic spectrum was obtained from a global fit to the flare data and 
consisted of two power-laws T(E S ) = (E s /lTeV)~ a with different spectral in- 
dices: a = 2.4 for the gaussian excess component and a = 2.7 for the contin- 
uous baseline [20]. These values were changed within few standard deviations 
without a significant impact on the final results. 

The energy resolution was assumed to be linear with energy <je = 0.22 x E 
although several other more complex energy dependencies consistent with the 
MAGIC telescope energy resolution measurements were tested with no signif- 
icant impact in the results. 
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Fig. 1. Parameterized flare shape for the Markarian 501 MAGIC flare data: the data 
are well described by a Gaussian-shaped flare activity on top of a constant baseline 
of steady gamma-ray emission (called "baseline" elsewere in this work). 



To infer from the data the best estimators for the fitting parameters, the 
likelihood function was fitted to the data using a standard numerical mini- 
mization package to minimize the function —2 log(L) as a function of four free 
parameters: 

• Mp/Mqcn the effective energy scale for a linear or a quadratic Quantum 
Gravity realization. 

• tw the gaussian intrinsic flare width. 

• to the maximum intrinsic flare position relative to the first gamma arrival 
time. 

• xb the normalization in area between the flat baseline component and the 
superimposed gaussian flare. 

The fit converged smoothly towards a minimum at Mqqi = O.SOjloao x 10 18 
GeV for the linear case and Mqg2 = 0.57loli9 x 10 11 GeV for the quadratic 
one. The uncertainties given above are asymmetric because they are given on 
Mg Gn and not on Mp/Mg Gn , the actual fitting parameter, and correspond to 
the change in the fitting parameter leading to an increase of A (—2 log(L)) = 1 
around the minimum. The shape of the Likelihood function as a function of 
Mp I Mqq\ around that minimum is almost parabolic as can be gleaned from 
Fig. [2] and therefore, the above uncertainty can be considered as enclosing 
one-sigma single parameter statistical uncertainty. Therefore, the fitted value 
of Mp/Mq G i differs from zero by over two-sigma. 
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Fig. 2. Behavior of the function —2 log(L) with respect to the parameter Mp/Mqgi 
around the minimum. The shape is rather parabolic and allows excluding the value 
Mp/Mqqi = at more than 2 sigma significance ( A(—2log(L)) = 4 ). 

The values of the fitting parameters in the linear and quadratic scan are shown 
in Tab. HJ As discussed in reference [29J this result is perfectly consistent 
with the one obtained by using the Energy Cost Function method described 
there. The correlation between the different parameters is always below 50% 
except between Mqgu and to, which are correlated at a level of 60% as the 
first parameter carries the information about the delay introduced by a QG 
effect during the propagation of photons and the second indicates where the 
maximum of the flare is located. 

The goal of this work is just the discussion of the new method described here 
and therefore we shall not discuss here the meaning or the interpretation of 
the above results. That is properly discussed in reference (29]. 

Moreover, a series of tests were performed to test the robustness, the stability 
and the interpretation of the result. 

First, the dataset was fitted to the same Likelihood function after scrambling 
randomly the arrival time of the photons. The fit then produced a best fit 
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n = 1 n = 2 



M P /M QGn 


o n+5.8 
°- y -4.5 


4.3l^ x 10 7 


tw 


219_2g s 


(228 ± 27) s 


to 


(2005 ± 42) s 


(2041 ± 36) s 


X B 


0.38 ± 0.04 


0.38 ± 0.04 



Table 1 

Best fit values of the free parameters in the Likelihood function for the linear and 
the quadratic scan. These values are obtained when all the parameters of the fit 
are left free, while the Mqg u values quoted in the text are obtained fixing all the 
parameters except for Mp/Mqgu itself. 

value compatible with zero energy-time correlation at one sigma. 

Second, mock simulated Monte Carlo datasets of the size of the actual data 
but with input parameters unknown by us were produced 2 and then fitted 
blindly by our Likelihood function, recovering the input parameters within the 
expected uncertainties. The same happened with the Energy Cost Function 
method as discussed in reference [29] . 

Third, thousand datasets of the same characteristics of the actual data were 
produced using a bootstrap method on the actual data. Each of these datasets 
was fitted to our Likelihood function and the distribution of the best-fit values 
was studied (see figure [3]). As expected from a probabilistic interpretation 
of the statistical uncertainties obtained in the Likelihood fit to the MAGIC 
data, the M QG1 = 0.30to; 2 Q x 10 18 GeV was enclosing about 68% of the best-fit 
values and less than 2% of the best-fit values were below zero. 

Therefore, it can be concluded that the Likelihood method explained here was 
producing in this reliable result. Moreover, reasonable modifications on 

the actual details of the physical process description casted in our Likelihood 
formula, as already discussed before, change the results in just a fraction of 
the statistical uncertainty. 

Finally, in order to check the goodness-of-fit of our Likelihood formula to the 
data, both the Likelihood function and the data have been binned following 
the binning originally used in reference [20], namely: four bins in energy and 
twelve bins in time (see Fig. H]). The x 2 obtained by comparing the data and 
fitting function histograms turns out to be in reasonable agreement with the 
number of degrees of freedom (equal to the number of used bins minus the 
number of fitting parameters). Therefore, as can be gleaned from Fig. HI our 
Likelihood function is indeed a good description of the MAGIC flare data. 

2 These datasets were independently produced by Adrian Biland with a toy Monte 
Carlo simulation implementing his own simple description of the physical process. 
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10 12 14 16 



M P /M QG1 



Fig. 3. Distribution for the best fit parameter Mp /Mqgx in bootstrap replica of the 
MAGIC flare data. As can be seen, just few fits produce Mp /Mqgi < 0. 

4 Example of application to simulated photons following a more 
complex flare structure 



The H.E.S.S. collaboration published in year 2007 the observation of a huge 
flare of the blazar PKS 2155-304 pi] , located at z = 0.116. The flare happened 
during the night of July 26th 2006, lasted about 90 minutes and was so intense 
(over 10 thousand gamma rays collected corresponding to a significance of 
almost 160 sigma) that a lightcurve with 1-minute time binning was produced. 

That lightcurve, shown in Fig. [51 exhibits a rather complex structure in time 
which, as discussed in reference [21], is well fitted by a set of five overlap- 
ping flares described by using "generalized Gaussian" shapes following the 
formulation suggested by Norris et al. [30J, that is: 

I(t) = Aexp[-(\t -t max \/a r4 ) K ] 

where a r (rise) is to be used for t < t max and cr^ (decay) is to be used for 
t > tmax- The actual best fit values for A, o r4 and k for each one of the five 



12 





Fig. 4. Comparison of the binned MAGIC data and Likelihood function distributions 
to check the goodness-of-fit. The overall % 2 /Ndf is 48.0/41. 
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Time - MJD53944.0 [min] 

Fig. 5. The lightcurve observed from the July 26th 2006 flare of PKS2155-304 by 
the H.E.S.S. collaboration (fig 1 in reference [H]). The data are binned in 1-minute 
intervals. The line is the fit to these data of the superposition of five bursts described 
in the text. 







Td 


K 


[min] [10~ 9 cm' 2 s" 1 ] 


N 


N 





41.0 


2.7±0.2 


173±28 


610±129 


1.07±0.20 


58.8 


2.1±0.9 


116±53 


178±146 


1.43±0.83 


71.3 


3.1±0.3 


404±219 


269±158 


1.59±0.42 


79.5 


2.0±0.8 


178±55 


657±268 


2.01±0.87 


88.3 


1.5±0.5 


67±44 


620±75 


2.44±0.41 



Table 2 

The results of the best x 2 fit °f the superposition of five bursts and a constant to 
the data shown in Figured! The constant term is 0.27 ± 0.03 x 10~ 9 cnr 2 s~ x (1.1 
Icrab)- Taken from [21]. 

overlapping flares are collected in Tab. El 

In addition, such a huge statistics allowed a precise analysis of the observed 
differential energy spectrum which, as can be seen in Fig. [6] turns out to be 
well described between about 200 GeV and 5 TeV by a broken power law such 
as: 



^ ^ dN 

E<EB] dE 
n „ dN 

e>Eb '1e 
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/o(T7^T7) r ^ ri ( 
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ITeV 
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Fig. 6. The time-averaged differential energy spectrum observed from the July 26th 
2006 flare of PKS 2155-304 by the H.E.S.S. collaboration (fig 3 in reference [21]). 
The dashed line is the best fit of a broken power-law to the data as described in the 
text. 

where J = (2.06 ± 0.16 ± 0.41) x 10~ 10 cm^s^TeV- 1 , E B = (430 ± 22 ± 80) 
GeV, T 1 = 2.71 ± 0.06 ± 0.10 and T 2 = 3.53 ± 0.05 ± 0.10 where for each 
parameter the two uncertainties are the statistical and the systematic values 
respectively. This energy spectrum has not been corrected for the absorption 
of VHE gammas on the Extragalactic Background Light but corresponds to 
the observed one, which is precisely the input needed for the correlation study. 

We wanted to test whether our method could be applied to such a kind of 
complex flare to eventually obtain results on any possible energy- dependent 
propagation delay effect. Nevertheless, since the method presented here relies 
on the use of the individual photon information and we have no access the 
actual H.E.S.S. data, Monte Carlo generated samples of simulated gamma- 
ray data following the H.E.S.S. published lightcurve and differential energy 
spectrum have been produced to try to understand the actual sensitivity our 
method would provide in spotting any non-trivial energy-time correlation, 
would it be present in the real data. 

For that, Monte Carlo simulated gamma-ray datasets with the same amount 
of gammas observed by H.E.S.S., their integral lightcurve, their differential 
energy spectrum and their claimed energy resolution were produced. In each 
dataset a different value for Mqc n was assumed: from no energy-time corre- 
lation (Mp/Mqcn = 0) to the same value found in the fits to the MAGIC 
data described before, and in the linear and in the quadratic scenarios. Those 
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datasets were fitted to a Likelihood function built following equation [T] in 
which the functions T(E a ) , G(E - E s ,a E (E s )) and F s (t - D{E„ M QGn , z)) 
are the same ones used for generating the simulated datasets, namely, the ones 
describing the H.E.S.S. data as discussed above. 

Fig. [7] shows the comparison between the binned Monte Carlo dataset and the 
Likelihood function for the case in which a QG scale equal to the one observed 
in the MAGIC data is assumed. In all tested cases, the fits do recover within 
one sigma the Mq Gn value used in the simulations, and this result demon- 
strates the reliability of our approach to analyze complex flare structures. 

Moreover, our fits predict that the H.E.S.S. data could provide a sensitivity 
in the determination of any measurement of Mqgu in the linear and in the 
quadratic assumptions, over a factor 6 better than the one obtained with the 
Markarian 501 flare observed by MAGIC. By playing with our fits, we've been 
able to trace back the origin of this improvement to three factors: 

• A factor about three comes from the fact that the redshift is about three 
times larger (0.116 compared with 0.034). 

• Another factor about two comes from the fact that the gamma statistics 
is almost ten times larger (about 12 thousand gammas compared with 1.4 
thousand gammas). That would produce a factor three smaller statistical 
uncertainty but, on the other hand, the observed MAGIC spectrum is harder 
(more photons at higher energies) than the H.E.S.S. one. In fact, it has 
been checked by playing with the spectral indices that the photons with the 
highest energies (above say 1 TeV) do play a crucial role in improving the 
sensitivity. 

• Finally an additional small gain comes from the rich flare structure although 
in general, the sensitivity is dominated by the fastest risetime or falltime in 
the whole lightcurve. 

We expect the H.E.S.S. collaboration to release soon the results of their anal- 
ysis of their flare and confirm or correct our expectations. 



5 Conclusions 

In this work we've presented a new approach for the analysis of correlations 
between the arrival time and energy of photons that can be applied to any 
observed set of photons coming from astrophysical sources to search for any 
kind of intrinsic or propagation energy-dependent delay effects. 

This approach uses directly the energy and arrival time information of each 
recorded photon and therefore does not require binning. Unlike other unbinned 
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Fig. 7. Comparison of the binned HESS Monte Carlo data and Likelihood function 
distributions to check the goodness-of-fit. The energy intervals chosen are the same 
ones used in fig. 2] whereas the time bins here are of 95 seconds duration. The overall 
X 2 /Ndf is 174.1/168. 17 



approaches, the one presented should be able to study the energy versus time 
correlations regardless on the complexity of the flare structure and even with 
quite low photons statistics. The method proposed is based upon a Likeli- 
hood estimator and therefore should be statistically optimal to find the best 
estimator for that correlation. 

This analysis technique was already applied to the Markarian 501 flare ob- 
served by the MAGIC telescope in July 2005 [201129] and some details about 
the actual implementation in that well as of several tests performed 

to verify the reliability, robustness and stability of the results have been given. 
For that case, the goodness-of-fit has also been discussed. 

Finally, the application to a more complex practical case such as the huge flare 
of PKS 2155-304 observed in July 2006 by the H.E.S.S. collaboration has been 
discussed. Monte Carlo simulations show that the method could be reliably 
used to analyze that data. 
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